% This function computes the specified centerized moments given the 
% quadrature weights over the continuous distributed grid points.
% INPUT:
%   QW_Mat(#continuous quadrature points x #discrete points):
%           Each column characterizes the quadrature,or histogram weight, 
%           over the continuously distributed dimensions, and each column
%           is corresponding to a non-continuously distributed dimension.
%           
%   QW_Unit (#grid points for continuous quadrature points x #continuous dimension):
%           Grid for continuously distributed quadrature points.
% 
%   OrderTable (#centerized moments x #dimension of continuous dimensions):
%           Table of centerized moment orders.
function [CentMom_Mat,Mean_Mat,Prob_Mat]=DistApp_QW_QW2Moment(QW_Mat,QQ_Unit,CentMomOrderTable)

%% Preliminary
% # of discrete dimension
N_Discrete  =   size(QW_Mat,2);
% # of continuous dimension
N_Continuous=   size(QQ_Unit,2);
% # of Moments
N_CentMom   =   size(CentMomOrderTable,1);

if size(CentMomOrderTable,2)~=N_Continuous
    error('Order Table does not match the continuous dimensions');
end

%% Compute and Collect the Moments
CentMom_Mat =   zeros(N_CentMom,N_Discrete);
Mean_Mat    =   zeros(N_Continuous,N_Discrete);
Prob_Mat    =   zeros(1,N_Discrete);

for ii=1:N_Discrete
    [Temp_CentMom,Temp_Mean,Temp_Prob] ...
                =   EmpiricalMoment(QQ_Unit,QW_Mat(:,ii),CentMomOrderTable);
    CentMom_Mat(:,ii) ...
                =   Temp_CentMom;
    Mean_Mat(:,ii) ...
                =   Temp_Mean;
    Prob_Mat(:,ii) ...
                =   Temp_Prob;
end

% MultiNum    =   size(QW_Mat,2);
% QQNum       =   size(QQ_Unit,2);
% MomNum      =   size(OrderTable,1);
% Moment_Mat  =   zeros(MomNum,MultiNum);
% Mean_Mat    =   zeros(QQNum,MultiNum);
% Prob_Mat    =   zeros(1,MultiNum);
% for ii=1:MultiNum
%     [CentMoment,Mean,TotalProb] ...
%             =   EmpiricalMoment(QQ_Unit,QW_Mat(:,ii),OrderTable);
%     Moment_Mat(:,ii) ...
%             =   CentMoment;
%     Mean_Mat(:,ii) ...
%             =   Mean';
%     Prob_Mat(:,ii) ...
%             =   TotalProb; 
% end
% 
% % Moment_Vec  =   reshape(Moment_Mat,[MomNum*MultiNum,1]);
% Prob_Vec    =   reshape(Prob_Mat,[MultiNum,1]);

return

function [CentMoment,Mean,TotalProb]=EmpiricalMoment(QQ,QW,OrderTable)

QQ_Dim          =   size(QQ,2);
QQ_Num          =   size(QQ,1);
Moment_Num      =   size(OrderTable,1);

TotalProb       =   sum(QW,1);
ConditionalQW   =   bsxfun(@rdivide,QW,TotalProb);

Mean            =   QQ'*ConditionalQW;

MaxOrder        =   max(max(OrderTable));
UnitCentMoment  =   zeros(QQ_Num,MaxOrder+1,QQ_Dim);
CentQQ          =   bsxfun(@minus,QQ,Mean');
for i=1:QQ_Dim
    UnitCentMoment(:,:,i)   =   bsxfun(@power,CentQQ(:,i),(0:MaxOrder));
end

CentMomentTable =   ones(QQ_Num,Moment_Num);
for i=1:Moment_Num
    for j=1:QQ_Dim
        Order_Temp      =   OrderTable(i,j)+1;
        CentMomentTable(:,i) ...
                        =   CentMomentTable(:,i).* ...
                            UnitCentMoment(:,Order_Temp,j);
    end
end

CentMoment      =   CentMomentTable'*ConditionalQW;



return